
STRATEGIES FOR ESTIMATING 

THE MARINE GEOID FROM ALTIMETER DATA 


P. Argentiero, W. D. Kahn, and R. Garza-Robles , 

Goddard Space Flight Center 
Greenhelt, Md. 20771 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. • 




TECH LIBRARY KAFB, NM 



TECH LIBRARY KAFB, NM 



□ 134 


S4 


1. Report No. 

NASA TN D-8285 

4. Title ond Subtitle 


2. Government Accession No. 


Strategies for Estimating the Marine Geoid from 
Altimeter Data 


7. Author(s) 

P. Argentiero, W. D. Kahn and R. Garza-Robles 

9. Performing Organization Name and Address 


Goddard Space Flight Center 
Greenbelt, Maryland 20771 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 


3. Recipient’s Catalog No. 

5. Report Date 

July 1976 

6. Performing Organization Code 

932 

8. Performing Organization Report No. 

G-7683 

10. Work Unit No. 

681-02-01-03 

11. Contract or Grant No. 


13. Type of Report and Period Covered 

Technical Note 

14. Sponsoring Agency Code 


15. Supplementary Notes 


16. Abstract 

In processing altimeter data from a spacecraft-borne altimeter to estimate the fine 
structure of the marine geoid, a problem is encountered. To describe the geoid fine 
structure, a large number of parameters must be employed, and to estimate all of 
them simultaneously is not possible. In practice, one is forced to hold a large 
number of parameters at a priori values and adjust others. Unless the parameteriza- 
tion exhibits good orthogonality in the data, serious aliasing results. Simulation 
studies show that, among several competing parameterizations, the mean free-air 
gravity anomaly model (i.e., Stokes’ formula) exhibited promising geoid recovery 
characteristics. 


Using covariance analysis techniques, this document provides quantitative measures 
of the orthogonality properties associated with the previously mentioned parameteri- 
zation. For instance, a 5- by 5-degree area mean free-air gravity anomaly can be 
estimated with an uncertainty of 0.01 mm/s 2 (1 mgal) (40-cm undulation) if all 
free-air gravity anomalies within a spherical radius of 10 arc degrees are simultane- 
ously estimated. 


|17. Kay Word* (S«l«ct«d by Author(s)) 


18. Distribution Statement 


Altimeter, Marine geoid. 
Gravity anomalies 


Unclassified— Unlimited 


|19. Security Classif. (of this report) 

Unclassified 


Cat. 46 


20. Security Classif. (of this page) 

21 . No. of Pages 

22. Price* 

Unclassified 

27 

$3.75 


For sale by the National Technical Information Service, Springfield, Virginia 22161 



This document makes use of international metric units according to the 
Systeme International d’Unites (SI). In certain cases, utility requires the 
retention of other systems of units in addition to the SI units. The conven- 
tional units stated in parentheses following the computed SI equivalents are 
the basis of the measurements and calculations reported. 


CONTENTS 


Page 

ABSTRACT i 

INTRODUCTION 1 

SKYLAB AND GEOS-C ALTIMETER EXPERIMENTS 2 

DESCRIPTION OF ALTIMETER MEASUREMENT 4 

ORTHOGONALITY AND ALIASING 15 

ESTIMATION STRATEGIES FOR GRAVITY ANOMALY 

DETERMINATION 20 

SUMMARY 24 

ACKNOWLEDGMENT 25 

REFERENCES 27 


iii 



STRATEGIES FOR ESTIMATING THE MARINE 
GEOID FROM ALTIMETER DATA 


P. Argentiero 
W. D.Kahn 
R. Garza-Robles 

Goddard Space Flight Center 
Greenbelt, Maryland 


INTRODUCTION 

The primary function of the spacecraft-borne altimeter is to determine the fine structure of 
the mean sea-surface topography. The instrument is well suited for the task. Consider, for 
instance, the altimeter on board the Geodetic Earth-Orbiting Satellite-C (GEOS-C) spacecraft 
which was launched in late 1974. In its global mode, the instrument is capable of producing 
a measurement every 4 kilometers along the subearth path of the satellite. This implies 
that, even assuming considerable data compression, it will be mathematically possible to 
extract topographic detail of 1 degree or less from such data. 

But practical difficulties must be overcome before the full potential of altimetry as a data 
type can be realized. To determine these difficulties, standard estimation techniques 
used to obtain sea-surface topography from altimetry data must be closely analyzed. Essen- 
tially, the problem is to reconstruct an analytic surface from direct but noisy observations 
of the surface taken at specified points. The obvious approach is to parameterize the surface 
in terms of coefficients which define a suitably dense set of functions in the space of two 
dimensional analytic functions and to recover the coefficients from the data by means of a 
standard minimum variance filter. The physical setting of the problem often suggests the 
proper parameterization. If not, an arbitrary family of analytic functions such as the set of 
two dimensional polynomials can be used. For this problem, a natural parameterization is 
suggested by the fact that, after suitable corrections, the mean sea surface can be considered 
as cohering closely to an equipotential surface known as the marine geoid (Reference 1). 
Thus, one candidate for a natural parameterization is the set of standard spherical harmonic 
coefficients of the Earth’s geopotential field. 

Any parameterization capable of describing the fine structure of the surface in question must 
employ an enormous number of coefficients. For instance, determinization of 3-degree 
features of the marine geoid using the standard spherical harmonic expansion of the geopo- 
tential field as a parameterization would require a full set of coefficients up to degree and 
order 60. This implies the estimation of over 3700 parameters. Unless special circumstances 
apply, it is impossible to simultaneously estimate such large parameter sets. In practice, 



to use parameter estimation techniques in recovering the marine-geoid fine structure from 
altimeter data, it will be necessary to adjust small subsets of the parameters while “freezing” 
the rest at a priori values. But, unless the parameterization exhibits a certain property with 
regard to altimetry data which is termed “orthogonality in the data set,” uncertainties in 
the adjusted terms will badly corrupt the estimates of the adjusted terms. This effect is fre- 
quently called aliasing. The orthogonality property just mentioned (defined in detail in a 
later section) is a property of a parameterization which permits a decomposition of the esti- 
mation problem into estimation problems of much smaller dimensionality without fear of 
serious aliasing. 

To take full advantage of the attractive properties of altimetry, it is necessary to develop a 
parameterization of the marine geoid which exhibits good orthogonality properties in altim- 
eter data. The orthogonality properties of the set of spherical harmonic coefficients of the 
geopotential field are poor and are not good candidates for the proper parameterization. 

Other candidates for the parameterization of the marine geoid are surface density layers 
(Reference 2) and sample functions (Reference 3). The parameterization whose properties 
are investigated in this paper is the one provided by mean free-air gravity anomalies and Stokes 
formula (References 4 and 5). If this parameterization is used to determine the marine geoid 
from altimeter data, it is important to determine to what extent mean free-air gravity 
anomalies are orthogonal in the data. Specifically, the important factor is how far two 
mean free-air gravity anomalies must be separated to ensure that neglecting one anomaly 
does not badly alias the estimate of the other. Without knowledge of this distance, intelligent 
use of the parameterization is impossible. 

Results of the Skylab altimeter experiment have demonstrated the ability of altimetry to re- 
veal marine-geoid fine structure. However, the GEOS-C mission will provide the first oppor- 
tunity for scientists to examine the output of a spacecraft-borne altimeter functioning in a 
global mode over the oceans of the world. The next section details the performance of the 
Skylab altimeter and the goals of the GEOS-C experiment. The mathematical details of 
the measurement modeling and preprocessing of altimeter data will be described in a later 
section, and the Stokes’ formula mean free-air gravity anomaly parameterization of the 
marine geoid will be documented. Following that will be a detailed discussion of the dual 
concepts of orthogonality and aliasing and their relationships to the problem of estimating 
a marine geoid from altimeter data. Finally, the results of a systematic application of co- 
variance analysis will be used to develop optimal estimation strategies for determining the 
marine geoid from altimetry data. 

SKYLAB AND GEOS-C ALTIMETER EXPERIMENTS 

The GEOS-C spacecraft was inserted into orbit during the latter part of calendar year 
1974. The prime experiment of this spacecraft was a radar altimeter. The objectives of 
the GEOS-C altimeter experiment were: (a) to determine the feasibility and utility of a 
spacebome altimeter to map the topography of the ocean surface with an absolute height 
accuracy of ±5 meters and with a relative height accuracy of 1 to 2 meters; (b) to 
determine the feasibility of measuring the deflections of the vertical at sea; (c) to determine 
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the feasibility of measuring wave height; and (d) to contribute to technology leading to a 
future operational satellite-altimeter system with a 10-cm measurement capability. 

When calibrated and corrected (e.g., for sea state, ocean tides, and other effects), the altimeter 
data constitute measures of the distance between the GEOS-C spacecraft and the ocean sur- 
face. Knowledge of the satellite altitude relative to a reference ellipsoid and knowledge of the 
oceanographic departures of the sea-surface topography from the geoid will then permit the 
determination of the geoid. The chief problem expected is the determination of orbital 
altitude for GEOS-C. The primary tracking systems for doing this are the satellite-to-satellite 
tracking system and precision lasers. Data from these and other systems tracking GEOS-C, 
including C-band radars, USB range and range-rate stations, and TRANET Doppler stations, 
will be used to find satellite height. 

Satellite contributions to the determination of the current ocean geoid have spatial resolutions 
corresponding to half wavelengths of approximately 18 degrees (2000 km). Although surface 
gravity achieves representations with finer resolution, in the 1- to 5-degree (1 10- to 550-km) 
range, it covers only about 50 percent of the ocean surface. The GEOS-C altimeter system 
will therefore fill in the gaps and provide valuable independent knowledge where data now 
exist. 

The Skylab altimeter has recently demonstrated the ability of a radar altimeter to detect 
features of the ocean surface. Figure 1 is an analysis of an altimeter pass over the Puerto 
Rican Trench. The pass was over the southwest corner of Puerto Rico and was 72.8 seconds 
in duration. A plot of the height residuals (formed by comparing the altimeter measure- 
ments with the calculated height of Skylab ’s orbit) shows that a 17-meter drop in height 
residuals occurred when Skylab passed over the deepest part of the ocean (i.e., corresponding 
to a 4000-meter depth of the ocean bottom) and that the peaking of the height residuals 
occurred when Skylab traversed the Puerto Rican land mass. The Skylab data described 
here exemplify the high level of resolution of surface features by a radar altimeter. 

Note, however, that data from the GEOS-C experiment will be significantly free of effects 
of spacecraft dynamic motions, (which is not the case for Skylab altimeter data because of 
attitude-control jet thrusting, crew motion, etc.). In addition, because GEOS-C had a nomi- 
nal 1-year operating lifetime, the altimetry data obtained was far greater than that obtained 
from Skylab. For reasons previously cited, the GEOS-C altimeter experiment data are more 
suitable for improving the marine geoid. 

To successfully use GEOS-C altimetry data for improving the accuracy of the marine geoid, 
the development of unique computer programs capable of processing these data was neces- 
sary. The following sections describe the mathematical models associated with these pro- 
grams and the computation strategies for processing altimeter data. 
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{BEGIN PASS) 




Figure 1. Altimeter pass over the Puerto Rico Trench: (a) ground track 
of altimeter pass; (b) Skylab-I altimeter residuals versus time data 
source: NASA/Wallops Station. 

DESCRIPTION OF ALTIMETER MEASUREMENT 

Figure 2 shows the geometry associated with the altimeter measurement. As can be seen, 
the altimeter measurement is nominally the shortest distance between the satellite and the 
sea surface (i.e., the measurement is along the normal to the sea surface that passes through 
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Figure 2. Altimeter measurement geometry. 


the satellite). The following relationship gives the mathematical model for the altimeter 
measurement: 

h a = h-N'-Sh - Ah' O) 

a u 

where 

h = ip-Tj 

p = position vector of S/C 

7 _ geocentric radius vector 

N' = geoid height above reference ellipsoid 

8h = deviation of sea surface from geoid 

S 

Ah = systematic errors in altimeter measurement (e.g., refraction, timing, etc.) 
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A more detailed description of the mathematical modeling of the altimeter measurement 
appears in Reference 6. 

The error sources that effect the altimeter measurement fall into three categories: (1 ) those 
that result from orbit uncertainty; (2) those that cause the measured geoid to deviate from 
the true geoid; and (3) those that affect the measurement accuracy itself. Equation 2 de- 
scribes the functional dependence of the error sources on the altimeter measurement: 

h a = h(E 0 t) - N'(G Q ) - 5h s (T o ,s o ,0,A,t) - Ah'(m 0 ,0,X,t) (2) 


where 

h = S/C altitude above reference ellipsoid 

E = orbit parameters and orbit dependent terms (radiation pressure, drag, etc.) 
t = time of altimeter measurement 
G 0 = vector of geopotential coefficients 
N' = geoid undulation function 
7 = vector of tidal error model coefficients 

1 = vector of local sea-surface biases (i.e., currents, storm surges, wind waves . . . ) 

m 0 = vector of altimeter measurement error coefficients 


or in general 


h a = h a (E o ,G o ,? o ,s o ,m o ,0,X,t) (3) 

The altimeter measurement will be corrected for known error sources, smoothed, sorted, 
etc., by the data preprocessor. Each measurement will then be compared with altitude cal- 
culated from the satellite’s orbit (h a ') to form the measurement’s residual Ah a 

Ah a = h a - h 'a (4) 

Mathematically, this residual is equivalent to the difference function, 
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or 


Ah a = Ah a (E) + Ah a (G) + Ah a (r) + Ah a (s) + Ah a (m) 


(5a) 


Models for the altimeter measurement errors follow. 

Error Attributable to Orbital Parameters (Ah a (E)) 

Precision laser tracking, satellite -to-satellite tracking (SST), USB, C-band, and Doppler data 
preceding and during altimeter measurement data acquisition will be used to correct the 
error on altimeter measurements attributable to orbital parameters using an existing orbit- 
determination program. Orbit-height accuracy required for altimetry data reduction must 
be better than 1 meter. Recent studies indicate that this is feasible (Reference 7). 

Error Attributable To Tides And Sea Surface (Ah g (r) and Ah a (S)) 

Corrections for deviations from the geoid attributable to ocean tides are based on the model 
of Hendershott (Reference 8). This model represents the dominant lunar semidiurnal tide 
(M 2 ) with lunar declination terms neglected and is representative of the current state of the 
art. The maximum error contribution attributable to tides is approximately 1 meter. The 
contributions attributable to local sea-surface biases are on the order of 50 cm, below the 
resolution threshold of the GEOS-C altimeter. Therefore, sea-surface bias corrections are 
neglected. 

Errors Attributable To Altimeter Hardware (Ah g (m)) 

Altimeter measurement errors which directly affect the accuracy for determining the geoid 
radius are: 

Ah a (m) = Ah 0 +Ahj + Ah L +Ah D +Ah x +Ah R (6) 

where 

Ah 0 = altimeter antenna offset 

Ahg = antenna offset due to S/C libration 

Ah L - dynamic lag 

Ah D = altimeter drift 

Ah T = timing bias 

Ah R = tropospheric refraction 

Models for the foregoing error sources follow. 
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Altimeter Antenna Offset 


Ah 0 = h a tan Vi(ol q - 5 b ) tan(a 0 - 5 b ) 
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where 

h a = spacecraft altitude 

5 b = antenna half beamwidth angle 

a 0 = angle off boresight 

r = transmitted pulse length 

c = velocity of light (299.7925 X 10 3 km/sec) 

Antenna Offset Attributable to Spacecraft Libration 

Ahg = h a tan V 2 (a% - 5 b ) tan(a £ - 5 b ) 


(cr) 3 

‘A 

[/A\ 2 - 22s] 

45 b h a_ 


lUJ J 


(6a) 


(6b) 


(6c) 
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«b < “s < 2° 
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(6d) 


for 

0 < a £ < 5 b 

where 



a 2 = S/C libration angle 

A e = peak libration angle 

</> e = libration phase angle 

T = orbital period 


Dynamic Lag 

Dynamic lag is a measurement-bias error induced by the change in spacecraft position at the 
time of the outgoing and return pulses. 


o o 



(6e) 


where 

O 

hj = altitude rate at time t. (sec) 
JC = lag coefficient (sec 2 ) 
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Altimeter Drift 

Altimeter drift is a measurement-bias error attributable to component aging. 

Ah D = D(t-t 0 ) (6f) 

where . 

t Q = initial time of altimeter measurement 
t = current time of altimeter measurement 
D = drift coefficient (sec) 

Timing Bias 

Timing bias is bias error introduced by spacecraft clock error. 

o 

Ah x = hr (6g) 

o 

h = height rate at time of altimeter measurement 
t = clock error coefficient 


Tropospheric Refraction Error 

Tropospheric refraction error is a bias error introduced into altimeter measurement because 
of tropospheric refraction. The model used is that developed by J. Saastamoinen (Reference 9). 

Error Attributable to Geopotential Ah a (G) 

In the preprocessing of altimeter data, the altimeter measurement will have been corrected for 
the orbit uncertainties of ocean tides, local sea-surface biases, and measurement-bias errors 
using the error models described previously. Thus equation 4 is restated as follows: 

Ah a (G) = h a - [h; + Ah a (E) + Ah a (r) + Ah a (s) + Ah a (m)] 

The relationship between altimeter measurement residuals to geoidal parameters can there- 
fore be stated as follows: 


Ah a (G) 



3G; 


AGj. 


( 8 ) 
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As stated previously, the altimeter measurement residual as described in Reference 7 
solely reflects primarily (neglecting second order effects) the departure of the geoid from 
the reference ellipsoid or, more precisely, the distance along the unit normal to the refer- 
ence ellipsoid between the reference ellipsoid and the geoid which is called geoidal height 
or geoidal undulation N' (figure 2). 

The mathematical relationship which relates the geoidal undulation to the disturbing poten- 
tial T is given by Bruns’ formula derived in Reference 4. That is, 

T 

N =— (9) 


T = disturbing potential 

7 = magnitude of gravity vector normal to reference ellipsoid 


The disturbing potential of the global geoid at point (0, X, r) is expressed in terms of spherical 
harmonics as follows: 


GM 'r'' /a\ n ^r— » 

r - — 2^ ( 7 ) 2^ P nm< sin *) l C nm c °s mX + s nm sin mX] (10) 


n= 2 


m = 0 


where 


GM 
C , S 

nm’ nm 
P nm ( §in *) 

(0, X) 


a 

. L 

• n 

m 


= gravitational constant of the earth 

= spherical harmonic expansion coefficients 

= associated Legendre functions 

~ latitude and longitude on geoid at which disturbing potential is evaluated 

= geocentric radius from earth center of mass to evaluation point of disturbing 
potential on geoid 

= semimajor azis of reference ellipsoid 

= is the limit of summation, and it is specified by the degree of harmonic 
expansion of the global geoid 

= summation index for degree terms of the spherical harmonic expansion 
of T 

= summation index for the order of terms in spherical harmonic expansion 
of T 
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Thus, the geoidal undulation at any point P (0, X, r) on the earth can be computed from geo- 
potential coefficients derived from satellites by analysis of perturbations on the orbits induced 
by the Earth’s gravity field. The undulations are computed from the combination of equa- 
tions 9 and 1 0. 

The Stokes’ formula (Reference 4, pp. 92-98) is another form of expressing the disturbing 
potential. This formula makes it possible to express the disturbing potential of the global 
geoid from gravity data. That is, 

T = 4®* / j AgS( ' l ' )da (U) 


where 


a 

= radius of reference ellipsoid 

f 

= flattening of reference ellipsoid 

S(« 

= Stokes’ function 

a 

= element of area 

Ag 

= surface gravity anomalies 


S(P) 



+ 1-6 cos ip - 3 cos 



$ , . 
— +sin 

2 


2 



(11a) 


from Bruns’ formula (equation 9). 

The geoidal undulation at any point P on the global geoid can be computed from Stokes’ 
formula (equation 1 1). That is' 


U- 


N = — =- — | ( Ag S(ip) da 

7 47T7 


(lib) 


In terms of geographical coordinates, Stokes’ function can be expressed as follows: 


N(0,A) = 


R 

47T7 


2ir /*7t/2 


" / \ ' _ n J j. ’ 


Ag(0',X') S(i p) cos (p'dtp'dX 1 


(11c) 


X' = 0 <p' = -ir/2 
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where 


da = cos 0'd0'dX' 

(0, X) = latitude and longitude of the computation point 
(0'X') = coordinates of the variable surface element a 

0 = spherical distance between the computation point and variable surface element 

0 = cos -1 [sin 0 sin 0' + cos 0 cos 0' cos (X - X')] 

Ag(0'X') = free-air gravity anomaly at the variable point (0'X') 

7 = mean value of gravity over Earth 

To combine surface gravity data and geopotential information derived from gravity field per- 
turbations acting on orbits of spacecraft, the Earth is divided into two areas (Reference 10), 
a global area (Aj) and a local area (A 2 ), surrounding the point P. Each gravity anomaly in 
each area is also partitioned into two parts represented by Ag s and Ag 2 respectively, where 

' L n n 

Ag s = 7 X X) (n ~ ') (7) P nm( sin ^{ C nm cosnX+S nm sin mX) (lid) 
_n = 2 m = 0 

The Ag 2 value is defined as the remainder of the gravity anomaly. By partitioning the Earth 
into two areas and the geoidal undulations into two components, equation 11c can there- 
fore be rewritten as follows: 

N(0,X) = N 1 +N 2 (lie) 

where 

^.2-n -rr/2 

N. = -5- I I Ag (0\X')S(0) cos 0'd0'dX' 

Jo J-7r/2 

or 

_ T 
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and 
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4ny 
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y: s {Ag s (0j Aj) } S( t//j ) COS 0j 
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(110 


where AN t is the correction to the geoidal undulations of the global geoid as a function of 
the corrections of mean free-air gravity anomalies. 

Equation 1 If is the form of the parameterization adapted for relating the altimeter measure- 
ment residuals to geoidal parameters. That is, equation 8 can now be written as follows: 


or in the matrix form 


5N = Ah = AN, 

« J 


5g = AGj = 5{Ag s (0jA])} 


A = 


3h R A^'AA' 


3 Gj 


4ttj 


S(^j) cos 0- 


5N (k X i) A (k X j) 5 §(j X l) 


( 12 ) 
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ORTHOGONALITY AND ALIASING 

Assuming that altimeter observations can be directly translated into deviations of the marine 
geoid from a reference geoid, it is straightforward to estimate gravity anomalies from altim- 
etry data. Let 5N be a vector of geoid undulations collected from a certain region over the 
ocean. Next, let Sg be a vector of mean free-air gravity anomalies defined over a region of 
the ocean which contains the data region. Then, repeating equation 12, 

5N = A5g (13) 

where A is a matrix the number of whose rows is equal to the number of data points and the 
number of whose columns is equal to the number of gravity anomalies. As demonstrated 
in equation 12, the individual elements of A (for example, A (I, J)) can be computed through 
Stokes’ formula and the latitude and longitude of the i* data point, together with the longi- 
tude and latitude of the midpoint of the grid over which the J* gravity anomaly is defined. 
Equation 13 provides a linear equation of condition, and, in a standard minimum variance 
fashion, 5g could be estimated from observations of 5N. However, in order for equation 13 
to be correct (correct, that is, assuming that the approximations inherent in the discrete 
form of Stokes’ formula are valid), the gravity anomalies in the array 5g must cover the 
globe. Considerably smaller regions would no doubt be adequate, but just how small these 
regions can be before a serious bias is introduced into the estimation process must be deter- 
mined by computations. In any case, computational considerations necessitate the choice 
of a region in which the number of elements in the Sg array does not exceed two or three 
hundred. Gravity anomalies outside this region are assumed to be zero. To see precisely 
what happens when this assumption is made, postulate that the Sg array of equation 13 is 
defined over an area so large that equation 13 is substantially correct. Then write 


5g = 



(14) 


where 

Sgj = gravity anomalies to be adjusted in a standard minimum variance filter 
and 

5g 2 = gravity anomalies assumed to be zero and thus left unadjusted by the minimum 

variance filter. 

Then equation 13 can be written 

6N = A t Sg'j + A 2 5g2 (15) 


15 


where Aj and A 2 are the variational matrices of 8 N with respect to 8 gj and with respect 
to 8 g^ , respectively. 

After proper corrections, altimeter data are treated as direct observations 8N of 8N with 
statistics 

8 N = SN + ^EOO - 0, E(w T ) = Q 

( 16 ) 

Estimates of mean free-air gravity anomalies obtained from satellite tracking and gravimetry 
measurements are available (Reference 10). Unless this information was correctly factored 
into the gravity anomaly estimates obtained from altimeter data, the resultant estimates 
would not be optimal. Consequently, we assume the existence of an a priori estimate 
8 gj ' of Sgj with properties 

8 g'j = 8 g L + oij , E (ttj) = U, E (a, a} ) = P t (n) 

For computational reasons, the gravity anomalies Sg 2 are assumed to be zero, but actual 
values of gravity anomalies in the region on the sphere which is ignored have a certain distri- 
bution about zero. We assume 


E( 8 g 2 ) = 0, E(Sg 2 Sgj) = P 2 


(18) 


When the values of Sg 2 are assumed to be zero, the minimum variance estimate of 6 g x be- 
comes 


Sg 1 = (Aj Q ' 1 Aj +P 7 1 )- 1 [Aj Q ” 1 8 N + P/ Sg' ] 


(19) 


See Reference 1 1 for details. Define the covariance matrix of the estimator given by equation 
19 as 

P = E [( 8 g 1 - 8 g 1 )( 8 g 1 -5gj) T ] (20) 

From equations 15, 16, 17, and 19, we obtain 

8 ^-Sgj = (Aj Q ' 1 A t +P - 1 )- 1 (-Aj Q ” 1 A 2 Sg 2 +Aj Q " 1 v + Pj 1 04 ) (21) 

Equation 2 1 yields 


P = (AjQ -1 Aj +P ^ l ) _1 + (Aj Q _1 Aj Pj 1 ) 1 Aj Q 1 A 2 P 2 AjQ 1 Aj 

(Aj Q _1 Aj +P 7 1 )- 1 


( 22 ) 
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Assume that the data is Uncorrelated and that each data point has the same variance. Then 

Q = I °l (23) 

where I is the identity matrix and o* is the common variance of the data. Also assume that 
the values of the unadjusted gravity anomalies are independently distributed. Then the co- 
variance matrix P 2 of 5g 2 Can be written as 


P 2 - 


0 


(24) 
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where n 2 is the number of unadjusted gravity anomalies and a? is the second moment about 
zero of the i* unadjusted gravity anomaly. 

Also define a matrix K as 


K = ( A J Q“ 1 A j +P 7 1 )- 1 A| CT 1 A 2 (25) 

If n t is the number of adjusted parameters, then K is of dimension n 2 by n 2 . With these 
assumptions, equation 22 yields the following expression for the variance of the i* adjusted 
gravity anomaly 

n 2 

P(U> - £ W,j of 

j = 0 (26) 

where /3? 0 is the i* diagonal element of matrix (A] - A ) _1 (assume here that diagonal elements 
of matrix P" 1 are relatively small) and 

j3 y = K(l,J), J > 1 (27) 

Define the error-sensitivity matrix as 

s = = 1.2,... ni ,J = 0,l,2,....n 2 (28) 

and, finally, define the alias matrix as 

L = so (29) 
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where 


a = 


o. 


(30) 


The alias matrix reveals much of the probability structure of the estimation procedure. 
Equations 26, 28, 29, and 30 show that the standard deviation of the i* adjusted gravity 
anomaly is the root sum square (RSS) of the terms in the i* row of the alias matrix. The 
elements in the first column of the alias matrix represent the RSS contribution to the stan- 
dard deviation of each estimated parameter caused by the data noise. The elements in the 
j* column, J > 2, represent the RSS contribution to the standard deviation of each esti- 
mated parameter attributable to the J - 1" unadjusted parameter. These terms are called 
the aliasing contributions to the uncertainty in the adjusted parameters attributable to the 
uncertainty in the J - l sf unadjusted parameter. Note that the aliasing contributions attri- 
butable to the J th unadjusted parameter are proportional to the standard deviation of the 
V h parameter. Definition: In a given estimation process, the \ th estim ated parameter is 
said to be orthogonal with respect to the J th unestimated parameter if the aliasing contri- 
bution to the i th estimated parameter attributable to the uncertainty of the i th unestimated 

parameter is zero. 

Two things can be noted about this definition. First, the orthogonality relationship between 
two parameters must be defined within the context of a given estimation process. Tire 
data set must be defined, and it must be stipulated which parameters are in the adjusted 
and the unadjusted modes. Second, although the discussion has been of mean free-air 
gravity anomalies and altimeter data, the mathematical results and definitions are applicable 
to any linear estimation problem in which some parameters are adjusted and others are 
left unadjusted. 

To see the implications of the orthogonality relationship, a more revealing representation of 
the aliasing terms is necessary. First, note that the first term on the right side of equation 
22 is the covariance matrix of the estimation process under the assumption that the un- 
adjusted parameters are perfectly known. This covariance matrix gives the uncertainty of 
the estimates attributable to data noise only. Define the so-called “noise-only” covariance 
matrix as 


P = (A}Q" 1 A 1 +P" 1 )" 1 


(31) 
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Next, observe that the elements in the i th row and 3 th column of A 1 and A 2 , respectively, 
are the partial derivative of the i th data point with respect to the 3 th adjusted parameter and 
the partial derivative of the i th data point with respect to the 3 th unadjusted parameter. The 
aliasing contribution to the i th adjusted parameter attributable to the 3 th unadjusted parameter 
can be written as 


n t m 

L(I,J+1)= £ P(I,K) 

K = 1 L = 1 

dSN(L) , 35N(L) 

3Sgj(K) ^ ( 95g 2 (J) 


(32) 


where m is the number of data points. 

If the estimates of the adjusted parameters are relatively uncorrelated in the noise-only 
covariance matrix, equation 32 can be approximated by 


m 

L(I,J + 1) = P(I,I) £ 
L = 1 


3SN(L) 95N(L) 

9§g j(I) ’ 95g 2 (J) 


(33) 


A sufficient condition for the left side of equation 2 1 to approximate zero is for the obser- 
vability patterns of 5gj (I) and Sg 2 (J) in the data to be virtually nonoverlapping. Up to 
about 40 degrees, Stokes’ function rapidly attenuates with increasing values of spherical 
radius (Reference 4). Hence, if the grids on which 5g, (I) and 5g 2 (J) are defined are suffi- 
ciently separated, the orthogonality relationship will be effectively satisfied and the estimate 
of Sg'j (I) will not exhibit aliasing from the uncertainty of 5g 2 (J). Conversely, if the grid 
on which 5gj (I) was defined were in close proximity to grids whose gravity anomalies were 
unadjusted, serious aliasing of the resultant estimate would be expected. 

It should be clear then, that if the gravity anomalies are estimated in a block, the outer layers 
of the block contain gravity anomalies whose estimates will be badly abased by the adjacent 
unadjusted parameters, and must therefore be discarded. However, the gravity anomalies in 
a sufficiently small inner core of the block may be adequately separated from the unadjusted 
parameters to be effectively orthogonal with respect to them. The estimates of these terms 
will presumably be of sufficient accuracy that they can be accepted. In effect, for every 
block of gravity anomahes that will be estimated, it will be necessary to construct a “buffer 
zone” several layers deep of gravity anomalies which surround the block. The new and larger 
block of gravity anomalies must be simultaneously estimated, and the estimates of gravity 
anomalies in the buffer zone must then be rejected because of aliasing. To achieve an 
intelligent compromise between estimation accuracy and computational load, it is necessary 
to determine the relationship between the depth of the buffer zone and the accuracy of the 
estimation procedure. The relationship will vary with grid size, computation radius, and 
data density. The most convenient way to study the relationship is to use covariance analysis 
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techniques to generate alias matrices for several situations and to attempt generalizations 
from the results. This is done in the following section. 

ESTIMATION STRATEGIES FOR GRAVITY ANOMALY DETERMINATION 

The accuracy with which a given gravity anomaly is estimated from altimeter data is a func- 
tion of many parameters. It is, of course, dependent on the accuracy and density of altim- 
eter data. As explained in the previous section, it is also dependent on the radial distance 
between the estimated gravity anomaly and the nearest gravity anomaly which is in the un- 
adjusted mode. This parameter is called the estimation radius. Another parameter, the 
computation radius, is an important determinant of the accuracy of a gravity anomaly esti- 
mation. The computation radius sets the maximum distance from a given grid element over 
which data is to be processed in order to estimate the gravity anomaly defined on that grid 
element. Estimation accuracy does not necessarily increase with increasing computation 
radius. Note that, for a given estimation radius, the covariance matrix of a set of estimated 
gravity anomalies is given by equation 22 as the sum of a matrix which is dependent only on 
data uncertainty and a matrix which represents the aliasing effects from the unadjusted 
parameters. With increasing computation radius, the elements of the first matrix must de- 
crease, but, in general, the elements of the matrix which conveys the aliasing effects will 
increase. This effect can be shown graphically by means of so-called aliasing maps. The 
aliasing maps of figures 3 and 4 were made simultaneously in a covariance mode, a situation 
in which the marine geoid was described by means of 2- by 2-degree gravity anomalies and 
by assuming 12 altimeter observations for each grid, each with an uncertainty of 1 meter. 

In figure 3, the RSS contribution was mapped in milligalileos (mgal)* to the uncertainty in the 
estimated gravity anomaly defined on the grid element in the lower left-hand corner when the 
adjacent anomalies were assumed to be in an unadjusted mode and to have an a priori uncer- 
tainty of 0.50 millimeters/sec 2 (mm/s 2 ) (50 mgal). The computation radius was 5 degrees. 
Note that the aliasing contributions decrease with the distance between the unadjusted 
parameter and the estimated parameter, demonstrating the inherent orthogonality property 
of the gravity anomaly parameterization of the marine geoid. In figure 4, the computation 
radius has been changed to 20 degrees, and the aliasing effect decreased much less radically. 

To achieve an intelligent compromise between computational load and estimation accuracy, 
it is necessary to determine the precise relationship between the accuracy of the estimate of 
a given gravity anomaly and the estimation radius and computation radius employed in the 
estimation procedure. To do this, it was assumed that altimeter data existed at a density of 
three data points per 1- by 1-degree grid and that the uncertainty on the data was 1 meter. 

It was also assumed that unadjusted gravity anomalies had a standard error about zero of 
0.50 mm/s 2 (50 mgal). This figure is conservative; a more realistic figure is 0.30 mm/s 2 
(30 mgal) (Reference 12). No a priori estimate was assumed for estimated parameters. The 
effect of unadjusted parameters out to a spherical radius of 45 degrees were included in the 


* Throughout the text of this document, the measurement unit of milligalileo (mgal) has been converted 
to a Standard International Unit of millimeters/sec 2 (mm/s 2 ). For convenience, the illustrations have 
not been converted. The simple conversion is 1 mgal = 0.01 mm/s . 
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computations. Under these conditions, figure 5 provides the standard deviation in the 
estimates of 5- by 5-degree mean free-air gravity anomalies as a function ot estimation 
radius and computation radius. This figure shows that the most efficient estimation strat- 
egy for 5- by 5-degree anomalies is obtained with an estimation radius of about 10 degrees 
and a computation radius of about 5 degrees. This strategy will result in an accuracy of 
0.01 mm/s 2 (1 mgal). This implies that if 5- by 5 -degree mean free-air gravity anomalies 
are estimated in blocks, estimates of gravity anomalies in the outer two layers of the block 
should be discarded because of aliasing. The rest of the estimated gravity anomalies should 
be accurate to within 0.01 mm/s 2 (1 mgal), provided a 5-degree computation radius is 
used. It is important to remember, however, that the prime intent is to estimate a marine 
geoid rather than gravity anomalies. With estimates of 5- by 5-degree mean free-air gravity 
anomalies, Stokes’ formula can be used to analytically reconstruct a marine geoid which is 
sufficiently detailed that 5 -degree features can be noted. Because the discrete form of 
Stokes’ formula is linear, it is easy to propagate a given error in gravity -anomaly estima- 
tion into an error in marine-geoid determination. Assuming that Stokes’ formula is accurate 
if its summation is carried out within a 90-degree spherical radius of a given point, an 0.01- 
mm/s 2 ( 1 -mgal) standard deviation in the estimate of 5- by 5-degree mean free-air gravity 
anomalies propagates into a 40-cm standard deviation of the resultant marine geoid. There- 
fore, by applying proper estimation procedures to the reduction of altimeter data, deter- 
mination of 5-degree features of the marine geoid with a resolution of 40 cm should be 
possible. 



Figure 5. Accuracy of 5- by 5-degree mean free-air gravity anomaly estimate versus compu- 
tation radius for various estimation radii. 
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To obtain a more detailed marine geoid presents a more difficult estimation problem. Figure 
6 provides the standard deviation in the estimates of 3- by 3-degree mean free-air gravity 
anomalies as a function of estimation radius and computation radius. An estimation radius 
of 12 degrees and a computation radius between 10 and 1 1 degrees appears to provide the 
most intelligent estimation strategy. The strategy leads to estimates of 3- by 3-degree mean 
free-air gravity anomalies with standard deviations of 0.05 mm/s 2 (5 mgal). This implies 
that, if 3-degree mean free-air gravity anomalies are estimated in blocks, the outer four 
layers must be discarded as being badly aliased. An 0.05-mm/s 2 (5-mgal) standard deviation 
in estimates of 3- by 3-degree mean free-air gravity anomalies propagates into a standard 
deviation of 1.2 m in the resultant marine geoid. With the present assumptions, no estima- 
tion strategy is adequate to determine features of the marine geoid as small as 2 degrees. 

With much greater data densities, such fine resolution is possible. For instance, if we in- 
crease the data density by a factor of 100, thus postulating 300 data points for each 1- by 
1-degree block, then, with an estimation radius of 10 degrees and a computation radius of 
2 degrees, 2- by 2-degree anomalies can be estimated with a standard deviation of 0.023 
mm/s 2 (2.3 mgal). This propagates into a standard deviation of the marine geoid of about 
42 cm. However, it is probably not realistic to postulate such data densities, and, in addi- 
tion, it suggests a very heavy computational load. The accurate resolution of 3-degree 
features of the marine geoid is likely to be the limit of what can be accomplished with 
altimeter data and the Stokes’ formula mean free-air gravity anomaly parameterization of 
the marine geoid. 



Figure 6. Accuracy of 3- by 3-degree mean free-air gravity anomaly estimate versus compu- 
tation radius for various estimation radii. 
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Finally, note that these results were accomplished without assuming a priori estimates of 
gravity anomalies. In fact, these estimates have been derived from satellite tracking data 
and direct gravity measurements. When these estimates are optimally combined with 
altimeter data by means of equation 10, the resultant estimates of gravity anomalies will be 
predicted on all relevant data. The representation of the marine geoid derived from these 
estimates will then constitute the best possible estimate of that surface. 

SUMMARY 

Data from the Skylab altimeter have demonstrated the ability of altimetry to discern the fine 
structure of the marine geoid. However, during the GEOS-C mission, global sensing of 
sea-surface topography from a spacecraft-borne altimeter will be possible for the first time. 

To make optimal use of the vast amount of altimeter data expected from the GEOS-C, 
numerous corrections must be made to the data, both to eliminate systematic errors inherent 
in the data type itself and to correct for the effects of deviations of the mean sea surface 
from the marine geoid. 

To employ standard parameter-estimation techniques for estimating a marine geoid from 
altimeter data, it is necessary to use a parameterization of the marine geoid which exhibits 
good orthogonality characteristics in the data type. In this respect at least, the Stokes’ formula 
mean free-air gravity anomaly parameterization appears to be the most promising. For a 
given grid size and data density, the accuracy of the geoid resulting from an estimate of 
gravity anomalies from altimeter data is a strong function of the choice of estimation radius 
and computation radius. An application of covariance analysis techniques reveals that, 
assuming a data density of three measurements per 1-by 1-degree spherical surface area with 
each measurement accurate to within 1 meter, optimal choices of estimation radius and 
computation radius of 10 and 5 degrees, respectively, yield an accuracy in the estimate of 
5- by 5-degree mean free-air gravity anomalies of approximately 0.01 mm/s 2 (1 mgal). This 
result implies that, with a judicious choice of estimation strategy, geoid height can be deter- 
mined with a standard deviation of 40 cm and a spatial resolution of 500 km (5 arc degrees). 
With the same data density and data accuracy assumptions, covariance analysis demonstrates 
that optimal choices of 12 degrees for estimation radius and 10 degrees for computation 
radius permits one to estimate with an accuracy of approximately 1 .2 m in geoid height with 
a spatial resolution of 300 km (3 arc degrees). Severe aliasing difficulties are encountered in 
attempting to estimate more detailed geoids. Our covariance analysis studies indicate, how- 
ever, that a very small computation radius, together with great data densities, can mitigate 
the aliasing effects and yield in localized areas a marine geoid capable of showing spatial 
resolutions at the 100- to 200-km level. 

Finally, note that use of the Stokes’ formula mean free-air gravity anomaly parameterization 
of the marine geoid makes it easy to combine satellite tracking data and direct measurements 
of gravity anomalies with altimeter data to obtain an optimal estimate of the marine geoid. 
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